#*************        丁香园-NHANES     *************
#*************        Analyze Code      *************
#*************       分析结果复现       *************
#*************                          *************

# Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
# Paper01-Dietary carotenoids and cognitive function among US adults, NHANES 2011–2014

#### 0.准备好环境 ####
library(gtsummary)
library(survey)
library(haven)
library(tableone)
library(plyr)
library(dplyr) # 链接：https://dplyr.tidyverse.org/reference/mutate.html
library(tidyverse)
library(arsenal) 



setwd("G:/BaiduNetdiskDownload/NHANES")
#setwd("J:/nhanes/数据分析/NHANES_20221011") #需要转换为自己的数据读取路径

#本文件是零代码论文适配文件-服务器配置使用文件

# 注：这里不需要引号，e.g. 直接写 CFDCSR 就可以了（视频中存在错误）-12.06更新
paper.data <- subset.data.frame(output, RIDAGEYR >= 60 &
                                  (!is.na(CFDCSR)) & #CERAD 延迟回忆评分
                                  (!is.na(CFDAST)) & #语言流畅性评分Animal Fluency
                                  (!is.na(CFDDS)) & #数字符号替代测试(DSST)
                                  (!is.na(DR1TLZ))& #饮食-叶黄素 & 玉米黄质摄取
                                  (!is.na(DR1TKCAL))& #能量摄入总量-Energy (kcal)
                                  (!is.na(RIDAGEYR)) & #年龄
                                  (!is.na(RIAGENDR)) & #性别
                                  (!is.na(RIDRETH3)) & #种族
                                  (!is.na(DMDEDUC2)) & #教育程度
                                  (!is.na(SMQ020)) #是否吸烟至少100支
)


dim(paper.data) #2712
colnames(paper.data)
analyze.variable <- c("SEQN", "WTDRD1", "SDMVPSU", "SDMVSTRA",
                      "Sex", "RIDAGEYR", "age.group", "race", "education.attainment", "INDFMPIR",
                      "alq.group", "smoke.group","BMXBMI", "BMI.group", "BMXWAIST", 
                      "total.calories", "total.dr.LZ", "total.ds.LZ", "total.dr.ds.LZ",
                      "CFDCST1", "CFDCST2", "CFDCST3",  "CERAD.total", 
                      "CFDCSR", "CFDAST", "CFDDS")


# 去除了分析的列
paper.data <- paper.data[, analyze.variable]
colnames(paper.data) <- c("SEQN", "WTDRD1", "SDMVPSU", "SDMVSTRA",
                          "Sex", 'Age', "Age.group", "Race", "Education.attainment", "PIR",
                          "Alq.group", "Smoke.group", "BMI", "BMI.group", "Waist",
                          "Total.calories", "Dietary.LZ", "Supplement.LZ", "Total.LZ",
                          "CERAD1", "CERAD2", "CERAD3", "CERAD.total",
                          "CERAD.delay.recall", "Animal.Fluency", "DSST")



#### 二、Paper Table ####
#### 1. 分析数据准备 ####
# 生成复杂抽样的对象
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust")



# 加权下的排序
# svytotal(~ Race, NHANES_design, na.rm=TRUE) # 离散性变量
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 


#### 2. Table1 右半截 #### 
tbl_svysummary(NHANES_design,  missing = 'no', 
               include = c(Age.group, Age, Sex,  Race, BMI.group, Alq.group, Smoke.group, Education.attainment,
                           Age, PIR, BMI, Waist, Total.calories, Dietary.LZ, Supplement.LZ, Total.LZ,
                           CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall, Animal.Fluency, DSST)) 

# 表格1 右侧非加权数据展示
tbl_svysummary(NHANES_design,  missing = 'no',
               include = c(Age.group, Sex,  Race, BMI.group, Alq.group, Smoke.group, Education.attainment,
                           Age, PIR, BMI, Waist, Total.calories, Dietary.LZ, Supplement.LZ, Total.LZ,
                           CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall, Animal.Fluency, DSST),
               statistic = list(all_continuous()  ~ "{median} ({p25}, {p75})", # 或者常用的："{mean} ({sd})"
                                all_categorical() ~ "{n_unweighted} ({p}%)"))%>% # 替换为非加权的 n
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)")



# 加权的情况，对连续变量进行离散分类
# https://rdrr.io/cran/survey/man/svyquantile.html
# 饮食中lz 的四分位数
Dietary.LZ.quantile.res <- svyquantile(~ Dietary.LZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1))
paper.data$Dietary.LZ.quantile.var <- cut(paper.data$Dietary.LZ,
                                          breaks = Dietary.LZ.quantile.res$Dietary.LZ[,'quantile'],
                                          labels = c('Q1', 'Q2', 'Q3', 'Q4'))

paper.data$Dietary.LZ.quantile.var[which(is.na(paper.data$Dietary.LZ.quantile.var))] <- 'Q1'

# 全部的lz 的四分位数
Total.LZ.quantile.res <- svyquantile(~ Total.LZ, NHANES_design, quantiles = c(0, 0.25, 0.5, 0.75, 1),na.rm = T)
paper.data$Total.LZ.quantile.var <- cut(paper.data$Total.LZ,
                                          breaks = Total.LZ.quantile.res$Total.LZ[,'quantile'],
                                          labels = c('Q1', 'Q2', 'Q3', 'Q4'))

paper.data$Total.LZ.quantile.var[which(is.na(paper.data$Total.LZ.quantile.var))] <- 'Q1'
#  # 将最小值替换为 Q1，以免为 NA
# 加权情况
# 注意：更新来 paper.data 后，需要更新 design 
NHANES_design <- svydesign(data = paper.data, ids = ~SDMVPSU, strata = ~SDMVSTRA, nest = TRUE, weights = ~WTDRD1, survey.lonely.psu = "adjust") 

##### 3.2 四分位数 Table2 一键生成结果 ####
# 使用函数 tbl_svysummary，来一键生成结果
# https://www.danieldsjoberg.com/gtsummary/reference/tbl_svysummary.html
# 表2 数据 第一行数据 Dietary L and Z intake
tbl_svysummary(NHANES_design, by = Dietary.LZ.quantile.var, missing = 'no',
               include = c(Dietary.LZ, Animal.Fluency, CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall,DSST),
               label = list(Dietary.LZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%               
  add_p()%>%
  as_flex_table() %>% 
  flextable::save_as_html(path = 'Table1_Result.html') 


# 显示表格2 下面 gcr Total L and Z intake
tbl_svysummary(NHANES_design, by = Total.LZ.quantile.var, missing = 'no',
               include = c(Total.LZ, Animal.Fluency, CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall,DSST),
               # label = list(Total.LZ ~ 'Total L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')
               ) %>%     # 设置变量的展示标签          
  modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%               
  add_p()%>%
  as_flex_table() %>% # 导出 html
  flextable::save_as_docx(path = 'Table2_Result.docx') 



# tbl_svysummary(NHANES_design, by = Dietary.LZ.quantile.var, missing = 'no',
#                include = c(Dietary.LZ, Animal.Fluency, CERAD1, CERAD2, CERAD3, CERAD.total, CERAD.delay.recall,DSST),
#                # label = list(Dietary.LZ ~ 'Dietary L and Z intake', Animal.Fluency ~ 'Animal Fluency: Score Total')
#                ) %>%     # 设置变量的展示标签
#   modify_header(all_stat_cols() ~ "**{level}**, N = {n_unweighted} ({style_percent(p)}%)") %>%
#   add_p() %>%
#   () %>% # 导出 html
#   flextable::save_as_html(path = 'Table3_Result.html')

